Introduction

Crucial to all good statistical analyses is a thorough exploratory analysis. In this document, we will introduce a few techniques for developing an understanding of our dataset, including an understanding of its limitations.

Load packages

Let’s load all the required R packages. See our tutorial 0_Installing_packages.html for installation information.

library(rgeos)
library(rgdal)
library(sp)
library(maptools)
library(spatstat)
library(spdep)
library(INLA)
library(inlabru)
library(readxl)
library(lubridate)
library(ggmap)
library(raster)

File set up

If you have not done so already, you need to download the two set of files associated with the workshop.

All the workshop script and tutorial files can be found on Github at github.com/joenomiddlename/DFO_SDM_Workshop_2020. To download these files on your computer, press on Code and then Download ZIP. Once downloaded, unzip the DFO_SDM_Workshop_2020-main.zip file. To be able to access some of the precompiled data, we need to make the unzipped folder the working directory. To do so in R Studio, navigate to the correct folder using the bottom right panel (folder view, you can use … button) and open it. Then, click “Set as Working Directory” under the tab ‘More’.

We should be inside the folder titled: ‘DFO_SDM_Workshop_2020’, you can verify this by using getwd().

Now we can load in the precompiled data

list2env(readRDS('./Data/Compiled_Data_new.rds'), globalenv())
## <environment: R_GlobalEnv>

Fin whale data

Throughout this workshop, we use data collected on fin whales in June, July, and August across the years 2007-2009 and for 2011. There are two main types of data:

  • sightings from systematic aerial surveys. These aerial surveys were conducted by both NOAA and DFO. Aerial survey sightings are accompanied with the aerial tracklines. Three surveys were combined: NOAA NARW Surveys (Cole, NOAA), 2007 T-NASS Aerial Survey (Lawson et al. 2009), and NOAA Cetacean Surveys (Cole et al., NOAA).

  • opportunistic sightings from two whale-watching operators. These sightings are maintained in the DFO Maritimes Region Whale Sightings Database (MacDonald et. al. 2017).

Coordinate reference system

These data are spatially-referenced data. The first thing we must always do is check for consistency between the coordinate reference systems (CRS) of each spatial object in use! To print the CRS of an sp object, simply add @proj4string to the name of the spatial object and run. For example, to check the CRS of the WW sightings we run:

Sightings_DRWW_sp@proj4string 
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs
We check the CRS of the other objects, but omit the code for brevity. 
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs

All of the spatial objects are in lat/lon - good! For future analysis we will be projecting the data into a different coordinate reference system to better preserve euclidean distance.

Finally, let’s turn off all warnings associated with coordinate reference systems. The latest PROJ6+/GDAL3+ updates have caused many warning messages to be printed.

rgdal::set_rgdal_show_exportToProj4_warnings(FALSE)
rgdal::set_thin_PROJ6_warnings(TRUE)
options("rgdal_show_exportToProj4_warnings"="none")

Plotting the data

Our first goal is generally to plot the data on a map. The gg() and gmap() functions from the inlabru package will prove extremely useful at plotting spatial data! The class of spatial objects we will use are from the sp package. The classes of these objects begin with ‘Spatial’. For example SpatialPointsDataFrame.

We have written a bespoke function gg.spatiallines_mod() to easily add SpatialLinesDatFrame objects to the plots too. This will prove useful for plotting transect lines. We load the bespoke functions to the working environment now.

source('utility_functions.R')

Let’s plot our data!

If the data are in lat/lon format then the gmap() function will automatically add a terrain layer from Stamen Maps to the plots. Let’s plot the survey sightings in blue, the survey tracklines in black, the whale watch sightings in purple, the whale watch ports in red.

gmap(Sightings_survey) +
  gg(Domain) +
  gg.spatiallines_mod(Effort_survey) +
  gg(Sightings_survey, colour='blue') +
  gg(Sightings_DRWW_sp, colour='purple') +
  gg(WW_ports, colour='red')

By default gmap uses Google maps (although this often changes (see ?get_map)), see Additional tips for ways to use open-source maps for publication.

To remove the map layer, simply replace the gmap(Sightings_survey) with ggplot().

ggplot() + # Notice the empty ggplot() call
  gg(Domain) +
  gg.spatiallines_mod(Effort_survey) +
  gg(Sightings_survey, colour='blue') +
  gg(Sightings_DRWW_sp, colour='purple') +
  gg(WW_ports, colour='red')

This plot hides some crucial information regarding the data collection. For example, the survey sightings and tracklines do not come from a single survey, or even a single organisation! Let’s plot this! The easiest way to do this is to subset the data accordingly!

table(Effort_survey$DATASET)
## 
##    DFO NOAA_1 NOAA_2 
##     49     54     70
# there are 3 surveys
ggplot() +
  gg(Domain) +
  gg.spatiallines_mod(Effort_survey[Effort_survey$DATASET=='DFO',], colour='purple') +
  gg.spatiallines_mod(Effort_survey[Effort_survey$DATASET=='NOAA_1',], colour='red') +
  gg.spatiallines_mod(Effort_survey[Effort_survey$DATASET=='NOAA_2',], colour='yellow')

This is problematic! The DFO tracklines (in purple) do not overlap with the two NOAA surveys! Thus, any future model will be unable to identify any differences in protocol efficiency. This is because any model intercepts will be confounded with the latent spatial field. More on this later!

Exercise 1

In addition, the surveys were conducted across 4 separate years. Let’s plot the survey tracklines by year. Do you see any cause for concern? Note that the names of the variables in the Effort_Survey are: DATASET and YEAR. Try this on your own! If you get stuck, click ‘Show Code’. Hint: The YEAR variable is of type character and contains 4 unique values (see below).

table(Effort_survey$YEAR)
## 
## 2007 2008 2009 2011 
##  101   20   19   33
class(Effort_survey$YEAR)
## [1] "character"
ggplot() +
  gg(Domain) +
  gg.spatiallines_mod(Effort_survey[Effort_survey$YEAR=='2007',], colour='purple') +
  gg.spatiallines_mod(Effort_survey[Effort_survey$YEAR=='2008',], colour='red') +
  gg.spatiallines_mod(Effort_survey[Effort_survey$YEAR=='2009',], colour='blue') +
  gg.spatiallines_mod(Effort_survey[Effort_survey$YEAR=='2011',], colour='yellow')
## Regions defined for each Polygons

Note that the surveys from years 2007, 2008 and 2009 covered largely different regions! Again, this is problematic if we want to model any changes in the whale distribution over time! The effect of year will be confounded by the spatial field. That being said, the data from 2011 appear to be a good candidate for model comparison as the spatial range overlaps with the other 3 years’ effort. We will holdout this data as our test data and use 2007, 2008, and 2009 as our training data.

xtabs(~ YEAR + DATASET, data=Effort_survey@data)
##       DATASET
## YEAR   DFO NOAA_1 NOAA_2
##   2007  49      3     49
##   2008   0     20      0
##   2009   0     19      0
##   2011   0     12     21

2011’s data comes exclusively from NOAA.

Transforming the data into a new CRS

For modelling, we will transform the data from lat/lon into a new “Canadian” CRS: NAD83 datum with UTM Zone 20N projection. This projection will help to preserve euclidean distance between points. The EPSG code is 2961.

To do the transformation, we will use the spTransform() function. For example, to transform the Whale Watch sightings spatial object Sightings_DRWW_sp, we first define the CRS object as follows:

Can_proj <- CRS("+init=EPSG:2961")
Can_proj <- fm_crs_set_lengthunit(Can_proj, unit='km')

The second line of code specifies that we want to work in units of km instead of the default meters. This can prove vital in applications to avoid numerical overflow.

Next, we transform Sightings_DRWW_sp:

Sightings_DRWW_sp <- spTransform(Sightings_DRWW_sp, Can_proj)
Sightings_DRWW_sp@proj4string
## CRS arguments:
##  +proj=tmerc +lat_0=0 +lon_0=-63 +k=0.9996 +x_0=500000 +y_0=0
## +ellps=GRS80 +units=km +no_defs
Notice the changed output from calling @proj4string. We repeat this for all the spatial objects that are points or lines. 
Sightings_survey <- spTransform(Sightings_survey, Can_proj)
Effort_survey <- spTransform(Effort_survey, Can_proj)
WW_ports <- spTransform(WW_ports, Can_proj)

Transforming the ‘raster’-like SpatialPixelsDataFrame objects (Bathym, Dist_Brier, and Dist_Quoddy) using spTransform would be inappropriate here. The projection leads to a curvature of the pixels. A more appropriate approach here is to use bilinear interpolation. The projectRaster() function from the raster package works great for this. This requires converting the SpatialPixelsDataFrame object into an object of type raster. This is made easy with the function raster(). Finally, to convert the raster object back into a SpatialPixelsDataFrame, we can use the as() function from the maptools package. This function is extremely useful for converting spatial objects between the popular packages: sp, spatstat, and sf. We use this function substantially throughout the workshop.

Bathym <- as(projectRaster(raster(Bathym), crs=Can_proj), 'SpatialPixelsDataFrame')
Domain <- spTransform(Domain, Can_proj)
Dist_Brier <- as(projectRaster(raster(Dist_Brier), crs=Can_proj), 'SpatialPixelsDataFrame')
Dist_Quoddy <- as(projectRaster(raster(Dist_Quoddy), crs=Can_proj), 'SpatialPixelsDataFrame')

Plot the (transformed) Bathymetry and Distance from Port spatial objects. We are going to combine these into a single plot using the multiplot() function from the inlabru package. This function takes as input ggplot objects and an argument layout, specifying how the plots should be arranged (see Additional tips for ways to change the layout).

multiplot(ggplot() +
  gg(Domain) +
  gg(Bathym) + xlab('East(km)') + ylab('North(km)') + labs(fill='Bathymetry') + 
    coord_fixed(ratio = 1),
ggplot() +
  gg(Domain) +
  gg(Dist_Brier) + xlab('East(km)') + ylab('North(km)') + 
  coord_fixed(ratio = 1),
ggplot() +
  gg(Domain) +
  gg(Dist_Quoddy) + xlab('East(km)') + ylab('North(km)')+ 
  coord_fixed(ratio = 1),
layout=matrix(1:4, nrow=2, ncol=2, byrow = TRUE))

Don’t like the colour scheme? See Additional tips to learn how to define your own manually.

Reducing the size of the domain to lie ‘close’ to the survey tracklines

In it’s current form, the domain is too large for the application. The domain extends far beyond the survey tracklines and the whale watch ports. To combat this, we will create a new domain that is restricted to lie within 30km of the nearest trackline.

# Create a set of SpatialPixels across the domain
pixels <- SpatialPoints(makegrid(Domain, n=5000),proj4string = Domain@proj4string)
# Convert the SpatialPixels into a SpatialPoints object before subsetting
#pred_int_points <- as(pred_int,'SpatialPoints')
# Subset to 30 km from Effort lines - can be slow
pixels_restricted <- 
  pixels[which(apply(gWithinDistance(pixels, Effort_survey,dist = 30,byid = T),2,sum)>0),]
# Select only the points that lie within the Domain
pixels_restricted <- pixels_restricted[Domain,]
# Convert to SpatialPolygons via SpatialPixels
pixels_restricted <- as(as(pixels_restricted,'SpatialPixels'),'SpatialPolygons')
# Merge the SpatialPolygons into a single Domain polygon using gUnionCascaded
pixels_restricted <- gUnionCascaded(pixels_restricted)
# Plot to check the new domain
ggplot() +
  gg(Bathym) + 
  gg(pixels_restricted) +
  xlab('East(km)') + ylab('North(km)') + labs(fill='Bathymetry') + 
  coord_fixed(ratio = 1)

# Smooth the Domain using gSimplify
Domain_restricted <- gSimplify(pixels_restricted,tol=20)
# Plot
ggplot() +
  gg(Domain_restricted) +
  gg(Domain, color='red') +
  gg.spatiallines_mod(Effort_survey, colour='blue')

Whilst the new SpatialPolygons object Domain_restricted no longer extrapolates far beyond the survey tracklines, the smoothing achieved by the function gSimplify() has caused the new Domain to no longer contain all the tracklines. To counter this, we will buffer the study region prior to simplifying using gBuffer.

Domain_restricted <- gSimplify(gBuffer(pixels_restricted,width=15),tol=20)
# Plot
ggplot() +
  gg(Domain_restricted) +
  gg(Domain, color='red') +
  gg.spatiallines_mod(Effort_survey, colour='blue')
## Regions defined for each Polygons

# Does the new domain contain all the survey tracklines?
gContains(Domain_restricted,Effort_survey)
## [1] TRUE

Great! We have defined a new buffered and smoothed domain that will no longer extrapolate far beyond our survey tracklines. To restore the original coastline definition we can use the function gIntersection(). We now update the SpatialPolygons object Domain using gIntersection():

Domain <- gIntersection(Domain_restricted,Domain)
ggplot() + gg(Domain) + gg.spatiallines_mod(Effort_survey)

Note that we will need to use Domain_restricted for defining our computational mesh later in the workshop! The smoothing of the coastline will help us to define a computational mesh that has fewer numerical issues. More on that later!

Extending the covariates to match the new domain

We will need to use Domain_restricted for defining our computational mesh later in the workshop! The package inlabru requires that our covariates are also defined at every point in the new domain. This will require extending our covariates to take ‘sensible’ values at points beyond our original domain (Domain). We will see that these imputed values will not affect the model estimates.

Our covariate depth is well-defined on the original (un-buffered) domain as defined by the SpatialPolygons Domain. At values outside this, we choose to ‘fill-in’ these points with nearest neighbour imputations. The code below does this for depth.

## Redefine log_Depth across our modified domain
# Create the log depth and log slope covariates
log_Depth <- Bathym
log_Depth$FIWH_MAR_Bathymetry[log_Depth$FIWH_MAR_Bathymetry >= 0] <- 0
log_Depth$FIWH_MAR_Bathymetry <- log(1-log_Depth$FIWH_MAR_Bathymetry)
names(log_Depth) <- 'log_Depth' 

# 1) Define a set of pixels across our modified domain
pixels_Domain <- as(SpatialPoints(makegrid(Domain_restricted, n=100000),proj4string = Domain@proj4string),'SpatialPixels')

# 2) Extract values of the covariate at the new pixel locations
pixels_Domain$log_Depth <- over(pixels_Domain,log_Depth)$log_Depth
# 3) impute missing values with the nearest neighbour value
pixels_Domain$log_Depth[is.na(pixels_Domain$log_Depth)] <- 
  log_Depth$log_Depth[nncross(as(SpatialPoints(pixels_Domain@coords[which(is.na(pixels_Domain$log_Depth)),]),'ppp'),
                          as(SpatialPoints(log_Depth@coords),'ppp'),
                          what = 'which')]
log_Depth <- pixels_Domain
# 4) Create a squared log depth covariate (for later)
log_Depth_sq <- log_Depth
names(log_Depth_sq) <- 'log_Depth_sq'
log_Depth_sq$log_Depth_sq <- log_Depth_sq$log_Depth_sq^2

# Plot the covariates with Domain_restricted overlayed
ggplot() + gg(log_Depth) + gg(Domain_restricted)

THE ABOVE HISTOGRAM SECTION WILL BE MOVED TO EXPLORATORY ANALYSIS

We assume the same functional relationship between each whale watch effort intensity \(\lambda_{eff, i}\) and distance from their port location \(s_i^*\). Here \(i\) indicates the port and \(s_i^*\) denotes the \(i^{th}\) port location. We scale the half-norm functions by a unique constant \(\lambda_i\), to capture the different number of vessels across the two ports: \[\lambda_{eff,i}(s) = \lambda_i exp\left(\frac{-||s-s_i^*||^2}{2\sigma_i^2} \right)\].

Dist_Surveylines$Dist_surveylines <- scale(Dist_Surveylines$Dist_surveylines^2)
Dist_Brier_sq <- as(Dist_Brier,'SpatialPixels')
Dist_Brier_sq$Dist_Brier_sq <- scale(Dist_Brier$Dist_Brier^2)
Dist_Quoddy_sq <- as(Dist_Quoddy,'SpatialPixels')
Dist_Quoddy_sq$Dist_Quoddy_sq <- scale(Dist_Quoddy$Dist_Quoddy^2)

Now we fit a Inhomogeneous Poisson point process model. As we are specifying the model on the log intensity, we need to include unique intercept parameters \(\alpha_{0,i} : i \in \{1,2\}\) for the two whale watch port intensities (\(\alpha_{0,i} = log \lambda_i\)) and three unique parameters \(\alpha_1\) and \(\alpha_{1,i} : i \in \{1,2\}\) defining the shapes of the three half norm functions (\(\alpha_{1,i} = \frac{-1}{2\sigma_i^2}\) and \(\alpha_1 = \frac{-1}{2\sigma^2}\)).

We must also inform spatstat about the source of each sighting, by fixing the value of its mark.

# All_obs_ppp$marks <- 
#   as.factor(c(rep('survey',dim(Sightings_survey@coords)[1]),
#             rep('Brier',dim(Sightings_Brier_nodup@coords)[1]),
#             rep('Quoddy',dim(Sightings_Quoddy_nodup@coords)[1])))

Ppm_2 <- ppm(All_obs_ppp~log_Bathym_im +
                    I(ifelse(marks=='Brier',1,0)):Dist_Brier_im +
                    I(ifelse(marks=='Quoddy',1,0)):Dist_Quoddy_im +
                    I(ifelse(marks=='survey',1,0)):Dist_Survey_im +
                    marks,
                  data=list(log_Bathym_im=log_Bathym_im,
                            Dist_Brier_im = as(as(Dist_Brier_sq,'SpatialGridDataFrame'),'im'),
                            Dist_Quoddy_im = as(as(Dist_Quoddy_sq,'SpatialGridDataFrame'),'im'),
                            Dist_Survey_im = as(as(Dist_Surveylines,'SpatialGridDataFrame'),'im'),
                            Domain=as(Domain, 'owin')),
                  subset = Domain)
coef(summary(Ppm_2))
##                                                       Estimate         S.E.
## (Intercept)                                       -977.9745269 123.64065586
## log_Bathym_im                                        0.4862812   0.08148644
## marksQuoddy                                        567.7490421 126.94313145
## markssurvey                                        697.5703053 133.08577596
## I(ifelse(marks == "Brier", 1, 0)):Dist_Brier_im   -819.1401759 104.13508711
## I(ifelse(marks == "Quoddy", 1, 0)):Dist_Quoddy_im -296.8180231  23.86004018
## I(ifelse(marks == "survey", 1, 0)):Dist_Survey_im -657.8434500 120.96619940
##                                                         CI95.lo      CI95.hi
## (Intercept)                                       -1220.3057594 -735.6432944
## log_Bathym_im                                         0.3265707    0.6459917
## marksQuoddy                                         318.9450764  816.5530078
## markssurvey                                         436.7269776  958.4136331
## I(ifelse(marks == "Brier", 1, 0)):Dist_Brier_im   -1023.2411962 -615.0391556
## I(ifelse(marks == "Quoddy", 1, 0)):Dist_Quoddy_im  -343.5828425 -250.0532037
## I(ifelse(marks == "survey", 1, 0)):Dist_Survey_im  -894.9328442 -420.7540558
##                                                   Ztest       Zval
## (Intercept)                                         ***  -7.909813
## log_Bathym_im                                       ***   5.967634
## marksQuoddy                                         ***   4.472468
## markssurvey                                         ***   5.241509
## I(ifelse(marks == "Brier", 1, 0)):Dist_Brier_im     ***  -7.866130
## I(ifelse(marks == "Quoddy", 1, 0)):Dist_Quoddy_im   *** -12.439963
## I(ifelse(marks == "survey", 1, 0)):Dist_Survey_im   ***  -5.438242

With this rough attempt at controlling for heterogeneous effort, we see that log depth is found to be significantly associated with the whale density, with a positive coefficient. The negative effects of all three complicated expressions imply that a negative association between ‘distance’ and observed intensity is witnessed (agreeing with common sense).

We can plot the estimated ‘true’ whale intensity \(\hat{\lambda}_{true}(s)\) predicted from this model:

# remove the effects of observers for plotting
pred_coef <- coef(Ppm_2)
pred_coef[c(1,3,4,5,6,7)] <- 0
plot(predict.ppm(Ppm_2, ngrid=c(300,300),
             new.coef=pred_coef)[[1]],superimpose=F,
     hcl.colors(60, "YlOrRd", rev = TRUE))

Alternatively, we can also plot the fitted observed intensity surface \(\lambda_{obs} =\lambda_{true} \lambda_{effort}\)) for all observer types. We plot the observed intensity surfaces for both the survey and Quoddy-based whale watch vessels.

pred_coef <- coef(Ppm_2)
plot(log(predict(Ppm_2, ngrid=c(300,300),
             new.coef=pred_coef)[[3]]),superimpose=F,
     hcl.colors(60, "YlOrRd", rev = TRUE))

# Plot Quoddy
plot(log(predict(Ppm_2, ngrid=c(300,300),
             new.coef=pred_coef)[[2]]),superimpose=F,
     hcl.colors(60, "YlOrRd", rev = TRUE))

We see that the observed intensity, \(\lambda_{obs}\) is dominated by the effort \(\lambda_{eff}\). Despite this, we still detect significant effects of both log depth.

This exploratory analysis is useful at informing us about possible species-environment relationships to investigate in our futurer models. By playing around with different covariates and functional forms (e.g. quadratic effects), we can develop an understanding of which models to formally compare in the modelling stage.

Note that the above effort-correction approach is not the same as the one performed by inlabru later. It should only be used for exploratory purposes.

Assessing the suitability of an inhomogeneous Poisson process

Unfortunately, by being in the inhomogeneous Poisson process class of point processes, the above models lack one crucial feature. They are unable to account for any additional spatial correlations/clustering that are frequently caused by unmeasured environmental covariates and biological processes acting behind the scenes.

This is why we turn to log-Gaussian Cox processes in the next tutorial. By assuming the intensity surface to be a realisation from a spatially-smooth (log-Gaussian) random field, we will be able to adjust our inference to take into account any additional spatial correlations and/or unmeasured covariates! By doing this, both parameter estimates and measures of uncertainty should be better adjusted to account for these autocorrelations. We will describe log-Gaussian Cox processes in depth in the next lecture.

Exercises

If you got stuck on any of the exercises, then please feel free to try them again. Here are links to the problems:

  1. Plotting data by year

Bonus Exercises

  1. Transform one of the spatial objects back into longitude latitude.
  2. Transform one of the spatial objects into units of meters.
  3. (More challenging) Investigate the relationship of distance from port with the two sources of whale watch data? What do you notice? What functional form could explain this type of behaviour? (Hint: use the over() function to extract the values of a spatial covariate (of type SpatialPixelsDataFrame) at point locations (stored as a SpatialPointsDataFrame object). Remember to subset the data accordingly!)

Additional tips and code

Using OpenStreetMaps instead of Google Maps

For publication, there can be issues regarding copyright of Google Maps. Using OpenStreetMap can help. To guarantee this simply add the following argument: source='osm', force=TRUE to gmap(). Double check the console that the maps are indeed being downloaded from stamen or osm. For brevity we have suppressed the messages.

gmap(Sightings_survey, source='osm', force=TRUE) +
  gg(Domain) +
  gg.spatiallines_mod(Effort_survey) +
  gg(Sightings_survey, colour='blue') +
  gg(Sightings_DRWW_sp, colour='purple') +
  gg(WW_ports, colour='red')

Note for this to work it needs to be in the original project (lat/lon), not UTM.

The multiplot() function is a very flexible function that enables publication-quality figures to be made with relative ease.

Changing the layout of multiplot()

To change the order of the plots you can change the argument byrow=TRUE to byrow=FALSE.

multiplot(ggplot() +
  gg(Domain) +
  gg(Bathym) + xlab('East(km)') + ylab('North(km)') + labs(fill='Bathymetry'),
ggplot() +
  gg(Domain) +
  gg(Dist_Brier) + xlab('East(km)') + ylab('North(km)'),
ggplot() +
  gg(Domain) +
  gg(Dist_Quoddy) + xlab('East(km)') + ylab('North(km)'),
layout=matrix(1:4, nrow=2, ncol=2, byrow = FALSE))

We can also change the size of the figures by changing the matrix in the layout.

multiplot(ggplot() +
  gg(Domain) +
  gg(Bathym) + xlab('East(km)') + ylab('North(km)') + labs(fill='Bathymetry'),
ggplot() +
  gg(Domain) +
  gg(Dist_Brier) + xlab('East(km)') + ylab('North(km)'),
ggplot() +
  gg(Domain) +
  gg(Dist_Quoddy) + xlab('East(km)') + ylab('North(km)'),
layout=matrix(c(1,1,2,3), nrow=2, ncol=2, byrow = TRUE))

As you can see plots of different size stretches the maps, to keep the set projection we can use coord_fixed(ratio = 1).

multiplot(ggplot() +
  gg(Domain) +
  gg(Bathym) + xlab('East(km)') + ylab('North(km)') + labs(fill='Bathymetry') +
    coord_fixed(ratio = 1),
ggplot() +
  gg(Domain) +
  gg(Dist_Brier) + xlab('East(km)') + ylab('North(km)') + 
  coord_fixed(ratio = 1),
ggplot() +
  gg(Domain) +
  gg(Dist_Quoddy) + xlab('East(km)') + ylab('North(km)') + 
  coord_fixed(ratio = 1),
layout=matrix(c(1,1,2,3), nrow=2, ncol=2, byrow = TRUE))

Define your own colour palette

We can define the colour palette easily.

colsc <- function(...) {
  scale_fill_gradientn(colours = rev(RColorBrewer::brewer.pal(11,"PuBuGn")),
                       limits = range(...))
}

Look at ?RColorBrewer::brewer.pal to see what other colour palettes are available.

multiplot(ggplot() +
  gg(Domain) +
  gg(Bathym) + xlab('East(km)') + ylab('North(km)') + labs(fill='Bathymetry') +
  colsc(Bathym@data[,1]),
ggplot() +
  gg(Domain) +
  gg(Dist_Brier) + xlab('East(km)') + ylab('North(km)') +
  colsc(Dist_Brier@data[,1]),
ggplot() +
  gg(Domain) +
  gg(Dist_Quoddy) + xlab('East(km)') + ylab('North(km)') +
  colsc(Dist_Quoddy@data[,1]),
layout=matrix(1:4, nrow=2, ncol=2, byrow = T))
## Warning in RColorBrewer::brewer.pal(11, "PuBuGn"): n too large, allowed maximum for palette PuBuGn is 9
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(11, "PuBuGn"): n too large, allowed maximum for palette PuBuGn is 9
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(11, "PuBuGn"): n too large, allowed maximum for palette PuBuGn is 9
## Returning the palette you asked for with that many colors

Have a go at creating your own colour palette function. Investigate the effects of changing both arguments to brewer.pal.

colsc2 <- function(...){
  scale_fill_gradientn(colours = rev(RColorBrewer::brewer.pal(7,"Spectral")),
                       limits = range(...))
}
multiplot(ggplot() +
  gg(Domain) +
  gg(Bathym) + xlab('East(km)') + ylab('North(km)') + labs(fill='Bathymetry') +
  colsc2(Bathym@data[,1]),
ggplot() +
  gg(Domain) +
  gg(Dist_Brier) + xlab('East(km)') + ylab('North(km)') +
  colsc2(Dist_Brier@data[,1]),
ggplot() +
  gg(Domain) +
  gg(Dist_Quoddy) + xlab('East(km)') + ylab('North(km)') +
  colsc2(Dist_Quoddy@data[,1]),
layout=matrix(1:4, nrow=2, ncol=2, byrow = T))

Potential repeated whale-watch and survey sightings in 2011

Could we not also put 2011’s whale-watch data into the training set? Below is some code showing that we cannot. A whale is sighted by at least one whale watch company on every day in 2011 that a survey detects a whale. The plot below shows that these sightings could have been of the same animal. Even if they weren’t, the tracklines in 2011 still visited areas that the whale watch vessels were present. Since our training and test datasets are required to be independent of each other, we must therefore also remove the 2011 whale watch sightings from the training data.

# what dates were survey sightings made on in 2011?
unique(Sightings_survey_test$DATE_LO[Sightings_survey_test$YEAR==2011])
## [1] "2011-06-21" "2011-08-12" "2011-08-13"
# Are sightings by the WW vessels made on these dates?
sum(grepl(Sightings_survey_test$DATE_LO[Sightings_survey_test$YEAR==2011][1],
     x=c(Sightings_Quoddy_nodup_test$WS_DATE[Sightings_Quoddy_nodup_test$YEAR==2011],
         Sightings_Brier_nodup_test$WS_DATE[Sightings_Brier_nodup_test$YEAR==2011])))>0
## [1] TRUE
sum(grepl(Sightings_survey_test$DATE_LO[Sightings_survey_test$YEAR==2011][2],
     x=c(Sightings_Quoddy_nodup_test$WS_DATE[Sightings_Quoddy_nodup_test$YEAR==2011],
         Sightings_Brier_nodup_test$WS_DATE[Sightings_Brier_nodup_test$YEAR==2011])))>0
## [1] TRUE
sum(grepl(Sightings_survey_test$DATE_LO[Sightings_survey_test$YEAR==2011][3],
     x=c(Sightings_Quoddy_nodup_test$WS_DATE[Sightings_Quoddy_nodup_test$YEAR==2011],
         Sightings_Brier_nodup_test$WS_DATE[Sightings_Brier_nodup_test$YEAR==2011])))>0
## [1] TRUE
# Could they be of the same animal? 
ggplot() + gg(Domain) + 
  gg(Sightings_survey_test[Sightings_survey_test$YEAR==2011,],colour='green') +
  gg.spatiallines_mod(Effort_survey_test[Effort_survey_test$YEAR==2011,],colour='yellow')

Acknowledgements

The code for hiding the Rmd code chunks came from Martin Schmelzer, found here

References to Data Sources:

References/Sources for data sets:

DFO Maritimes Region Whale Sightings Database

MacDonald, D., Emery, P., Themelis, D., Smedbol, R.K., Harris, L.E., and McCurdy, Q. 2017. Marine mammal and pelagic animal sightings (Whalesightings) database: a users guide. Can. Tech. Rep. Fish. Aquat. Sci. 3244: v + 44 p.

NOAA NARW Surveys

Timothy V.N. Cole. National Oceanic Atmospheric Administration, National Marine Fisheries Service, Northeast Fisheries Science Center. 166 Water Street, Woods Hole, MA, USA

2007 TNASS DFO Aerial Survey

Lawson. J.W., and Gosselin, J.-F. 2009. Distribution and preliminary abundance estimates for cetaceans seen during Canada’s marine megafauna survey - A component of the 2007 TNASS. DFO Can. Sci. Advis. Sec. Res. Doc. 2009/031. vi + 28 p.

NOAA Cetacean Surveys

Timothy V.N. Cole and D. Palka. National Oceanic Atmospheric Administration, National Marine Fisheries Service, Northeast Fisheries Science Center. 166 Water Street, Woods Hole, MA, USA

  1. Alternatively, we can use the function over() from the sp package. We choose to create im objects so that we can later use the ppm function from spatstat for exploratory analysis purposes. ppm fits Poisson process models, which can be equivalent to Maxent models.↩︎